function SVMRun
CreateGUI;
end

function CreateGUI(varargin)
%% Start GUI
persistent gui

diary(['diary/' num2str(datenum(clock)) datestr(clock,1) '.txt']);

gui.nobuttons        = [7 4];
gui.buttondistance   = [20 10];
gui.buttonsize       = [120 45];
gui.figuresize       = (gui.buttonsize+gui.buttondistance).*(gui.nobuttons)+gui.buttondistance;
gui.uipositions(:,1) = sort(repmat(round(linspace(gui.buttondistance(1),gui.figuresize(1)-gui.buttondistance(1)-gui.buttonsize(1),gui.nobuttons(1))'),gui.nobuttons(2),1));
gui.uipositions(:,2) = repmat(round(linspace(gui.buttondistance(2),gui.figuresize(2)-gui.buttondistance(2)-gui.buttonsize(2),gui.nobuttons(2))'),gui.nobuttons(1),1);

if ~isempty(varargin) && isa(varargin{1},'double')
    gui.source = varargin{1}; varargin={};
end
if ~isfield(gui,'fid') || ~ishandle(gui.fid);
    gui.fid          = figure(1);
    set(gui.fid,'position',[50 50 gui.figuresize]); else figure(gui.fid);
end

% Builds the gui controls
if ~isfield(gui,'source') || ~ishandle(gui.fid);
    gui.h.experiment        = uicontrol('String','SVM Gijs|SVM Bham|Test','Value',1,'Style','popupmenu','Callback',{@CreateGUI});
    gui.h.subjects          = uicontrol('String','Placeholder','Style','popupmenu','Callback',{@CreateGUI});
    gui.h.SVMtype           = uicontrol('String','winSVM|Matlab|Weka|Mex','Style','popupmenu','Callback',{@CreateGUI});
    gui.h.kerneltype        = uicontrol('String','polynomial','Style','popupmenu');
    gui.h.analysistype      = uicontrol('String','Placeholder','Style','popupmenu','Callback',{@CreateGUI});
    gui.h.shuffletimecourses= uicontrol('String','yes|no|both|','Value',3,'Style','popupmenu'); % Runs the program twice, with and without machine learner
    
    gui.h.plottimecourses   = uicontrol('Value',0,'Style','checkbox');
    gui.h.plotSVMresults    = uicontrol('Value',1,'Style','checkbox');
    gui.h.randomdata        = uicontrol('Value',0,'Style','checkbox'); % Put random data into SVM
    gui.h.parallelcomputing = uicontrol('Value',0,'Style','checkbox'); % Allow parallel computing (requires parallel toolbox)
    gui.h.exportmatrix      = uicontrol('Value',0,'Style','checkbox');
    
    gui.h.timecourses       = uicontrol('String','Temp','Style','popupmenu','Callback',{@CreateGUI}); % Shows available timecourses, select timecourse using timecoursestoanalyse (i.e. 1 = first tc in list)
    gui.h.smoothing         = uicontrol('String',num2str(10),'Style','edit');  % Smooth using moving avg. windowsize in seconds, Gijs used 10
    gui.h.bandpass          = uicontrol('String',num2str([0 20]),'Style','edit');  % Pass frequencies between [min max]
    gui.h.RespDelayCor      = uicontrol('String',num2str(200),'Style','edit'); % Response delay
    gui.h.normalization     = uicontrol('String','mean block, std scan|mean scan, std scan|mean','Style','popupmenu');
    gui.h.TimeShift         = uicontrol('String',num2str(0),'Style','edit'); % Shift the entire session in time (e.g. if fMR's are misalligned), in ms
    gui.h.CombineEvents     = uicontrol('String',num2str([-1 9]),'Style','edit'); % Amount of event combined in single volume after buttonpress (Event related)
    gui.h.NoLogs            = uicontrol('String',num2str(12),'Style','edit'); % Number of logs present
    gui.h.Tr                = uicontrol('String',num2str(2000),'Style','edit'); % Tr of fMRI scan in ms
    gui.h.onset             = uicontrol('String',num2str(20),'Style','edit'); % Length [sec] of extra period before start of trial (to avoid onset effect during data.hrf convolution
    gui.h.firstvol          = uicontrol('String',num2str(0),'Style','edit'); % First usefull fmr volume of every session (drop at least until this volume)
    gui.h.fmrvolumestoskip  = uicontrol('String',num2str(0),'Style','edit');
    gui.h.SVMruns           = uicontrol('String',num2str(12),'Style','edit'); % Number of SVM passes, must be even number!
    gui.h.trainsetsize      = uicontrol('String',num2str(4/8),'Style','edit'); % Size of training set (fraction)
    
    gui.h.start             = uicontrol('Style','pushbutton','String','Go','Callback',{@configuration});
end

%% Set dynamic fields

switch get(gui.h.experiment,'value')
    case 1; gui.subjects         = {'All','LD','TA','MJ','MM','TJ','TJFCD','TK','Specify'};
        gui.analysistypes    = {'FMr matrix','Event related','Multiclass','Beta','BetaAVG' ,'BetaState'};
        set(gui.h.trainsetsize,'String',num2str( 0.9));
        
        set(gui.h.SVMruns,'String',num2str(12));
        set(gui.h.NoLogs,'String',num2str(12));
    case 2; gui.subjects         = {'All','AM','HB','JZ','MJ','MP','MR','RW','Specify'};
        gui.analysistypes    = {'Blockedfewfeatures','Blocked'};
        set(gui.h.trainsetsize,'String',num2str( 0.75));
        
        set(gui.h.SVMruns,'String',num2str(12));
        set(gui.h.NoLogs,'String',num2str(4));
    case 3; gui.subjects         = {'TB','MJ','TG'};
        gui.analysistypes    = {'Test','Blocked','Blockedfewfeatures','FMr matrix','Event related','FMr matrix'};
        set(gui.h.SVMruns,'String',num2str(1));
        set(gui.h.trainsetsize,'String',num2str( 0.5))
        set(gui.h.NoLogs,'String',num2str(4));
end

set(gui.h.subjects,'String',gui.subjects);
set(gui.h.analysistype ,'String',gui.analysistypes);

switch get(gui.h.SVMtype,'value'); % winSVM|Matab|Weka|Mex
    case 1; allowedkernels = {'linear','dot','radial','neural','anova'};
    case 2; allowedkernels = {'polynomial','radial','linear'};
    case 3; allowedkernels = {'polynomial'};
    case 4; allowedkernels = {'polynomial'};
end;
set(gui.h.kerneltype,'String',allowedkernels)

analysistypes = get(gui.h.analysistype,'string');
switch char(analysistypes(get(gui.h.analysistype,'value')))
    case 'FMr matrix';        cfg.timecourses  = {'All','hrfbinStateTimecourse','hrfbinWantedTimecourse','hrfbinHoldSwitchTimecourse','hrfbinFlipTimecourse','hrfbinWantedflipTimecourse','hrfbinUnwantedflipTimecourse','hrfbinFlipToCCWTimecourse','hrfbinFlipToCWTimecourse','Specify'};
    case 'Multiclass';        cfg.timecourses  = {'All','LLvsLRTimecourse','RLvsRRTimecourse','LLvsRLTimecourse','LRvsRRTimecourse','LLvsRRTimecourse','LRvsRLTimecourse'};
    case 'Event related';     cfg.timecourses  = {'All','FlipERTimecourse','FlipERWantedTimecourse','FlipERUnwantedTimecourse','FlipERToCWTimecourse','FlipERToCCWTimecourse','FlipERWillTimecourse','FlipERWillToCWTimecourse','FlipERWillToCCWTimecourse','Specify'};
    case 'Blocked';           cfg.timecourses  = {'All','Memorysigned','Perceptionsigned','Halves','Specify'};
    case 'Blockedfewfeatures';cfg.timecourses  = {'All','Memorysigned2','Perceptionsigned2','Halves2','Specify'};
    case 'Test';              cfg.timecourses  = {'Test'};
    case 'Beta';              cfg.timecourses  = {'All','WantedvsUnwantedTimecourse','WillCWvsCCWTimecourse','Specify'};  
    case 'BetaAVG';           cfg.timecourses  = {'All','hrfbinStateTimecourse','hrfbinWantedTimecourse','hrfbinHoldSwitchTimecourse','hrfbinFlipTimecourse','hrfbinWantedflipTimecourse','hrfbinUnwantedflipTimecourse','hrfbinFlipToCCWTimecourse','hrfbinFlipToCWTimecourse','Specify'};
    case 'BetaState';         cfg.timecourses  = {'All','hrfbinStateTimecourse','hrfbinWantedTimecourse','hrfbinHoldSwitchTimecourse','hrfbinFlipTimecourse','hrfbinWantedflipTimecourse','hrfbinUnwantedflipTimecourse','hrfbinFlipToCCWTimecourse','hrfbinFlipToCWTimecourse','Specify'};
    
end
set(gui.h.timecourses,'String',cfg.timecourses)

% Modifies GUI positions
gui.uifields          = fieldnames(gui.h);
gui.checkboxcount     = 1; n=1; m=1;
while n<=length(gui.uifields);
    switch get(gui.h.(char(gui.uifields(n,:))),'Style')
        case 'checkbox'
            set(gui.h.(char(gui.uifields(n,:))),'String',char(gui.uifields(n,:)),'position',[gui.uipositions(m,1) gui.uipositions(m,2)+gui.buttonsize(2)*mod(gui.checkboxcount+1,2)/2 gui.buttonsize(1) gui.buttonsize(2)/2]);
            m=m-mod(gui.checkboxcount,2); gui.checkboxcount = gui.checkboxcount+1;
            if ~isequal(get(gui.h.(char(gui.uifields(n+1,:))),'Style'),'checkbox')&&mod(m,2)~=0; m=m+1; end
        case {'edit'}
            set(gui.h.(char(gui.uifields(n,:))),'position',[gui.uipositions(m,:) gui.buttonsize(1),gui.buttonsize(2)/2]);
            uicontrol('Style','text','String',char(gui.uifields(n,:)),'position',[gui.uipositions(m,1) gui.uipositions(m,2)+gui.buttonsize(2)/2 gui.buttonsize(1) gui.buttonsize(2)/2],'Backgroundcolor',get(gui.fid,'color'))
        case {'popupmenu'}
            set(gui.h.(char(gui.uifields(n,:))),'position',[gui.uipositions(m,:) gui.buttonsize(1),gui.buttonsize(2)/2]);
            uicontrol('Style','text','String',char(gui.uifields(n,:)),'position',[gui.uipositions(m,1) gui.uipositions(m,2)+gui.buttonsize(2)/2 gui.buttonsize(1) gui.buttonsize(2)/2],'Backgroundcolor',get(gui.fid,'color'))
            curstrings = get(gui.h.(char(gui.uifields(n,:))),'String');
            if isequal('Specify',char(curstrings(get(gui.h.(char(gui.uifields(n,:))),'Value'),:))) && isfield(gui,'source') && gui.source==gui.h.(char(gui.uifields(n,:)));
                gui.s.(char(gui.uifields(n,:))) = specifygui(gui.h.(char(gui.uifields(n,:)))); end
        case 'pushbutton'
            set(gui.h.(char(gui.uifields(n,:))),'position',[gui.uipositions(m,:) gui.buttonsize]);
        otherwise
            set(gui.h.(char(gui.uifields(n,:))),'position',[gui.uipositions(m,:) gui.buttonsize]);
    end
    set(gui.h.(char(gui.uifields(n,:))),'Backgroundcolor',get(gui.fid,'color'));
    n=n+1; m=m+1;
end

% Checks for nogui toggle
if ~isempty(varargin); for n=1:2:length(varargin)
        if ~isequal(varargin{n},'nogui')
            set(gui.h.(char(varargin{n})),'Value',varargin{n+1}); end; end;
end

gui.nogui = false;
configuration(1,{gui,cfg});

if sum(strcmp('noGUI',varargin))==1;         % GUI off
    gui.nogui = true;
    set(gui.fid,'Visible','off');
    configuration(gui.h.start,1);
end
end

function configuration(source,event)
%% Reads out GUI/sets variables

persistent gui
persistent cfg

if ispref('FMRISVM','datadir')
    cfg.datadir = getpref('FMRISVM','datadir');
else
    cfg.datadir = '';
end

cfg.LLcol=[0 1 0];
cfg.RLcol=[1 0 0];
cfg.LRcol=[0 1 1];
cfg.RRcol=[1 0 1];

disp(['datadir: ' cfg.datadir]);
 
if source == 1;
    gui                      =  event{1};
    cfg                      =  event{2};
    gui.strings.numberofruns = 0;
else
    
    cfg.CCWcode = -1;           % Counterclockwise rotation is indicated by this number. CW = -1*CCW.
    cfg.Holdcode = 1;           % Holding is indicated by this number. Switch =-Hold
    cfg.timestarted     = datestr(clock,1);
    cfg.uifields        = fieldnames(gui.h);
    
    for n=1:length(gui.uifields);
        switch get(gui.h.(char(gui.uifields(n,:))),'Style')
            case 'checkbox'
                cfg.(char(gui.uifields(n,:))) = get(gui.h.(char(gui.uifields(n,:))),'Value');
            case 'edit'
                z=1; remain = {};      remain = get(gui.h.(char(gui.uifields(n,:))),'String');
                while true;
                    [tefield remain]              = strtok(remain);
                    cfg.(char(gui.uifields(n,:)))(z) = str2num(tefield);
                    if isempty(remain); break; end; z=z+1; end;
            case 'popupmenu'
                
                structfields = cellstr(get(gui.h.(char(cfg.uifields(n,:))),'String')); value=get(gui.h.(char(gui.uifields(n,:))),'Value');
                if strcmp('Specify',structfields(value))
                    cfg.(char(gui.uifields(n,:))) = gui.s.(char(gui.uifields(n,:)));
                elseif strcmp('All',structfields(value))
                    cfg.(char(cfg.uifields(n,:))) = structfields(~strcmp('All',structfields)&~strcmp('Specify',structfields));
                else
                    cfg.(char(cfg.uifields(n,:))) = structfields(get(gui.h.(char(cfg.uifields(n,:))),'Value'));
                end
                
            case 'pushbutton'
                cfg.(char(gui.uifields(n,:))) = get(gui.h.(char(gui.uifields(n,:))),'Value');
            otherwise
                cfg.(char(gui.uifields(n,:))) = get(gui.h.(char(gui.uifields(n,:))),'Value');
        end
    end
    
    Runit(cfg,gui)
end
end

function Runit(cfg,gui)
%% Main body
cfg.NoOfEmptyLines = 8;     % Amount of header lines in log file (Check logfile syntax before running!
cfg.SampleRate = 100;       % M-seconds per sample (gets altered 1 sample to per fMRI volume later on)


for m = 1:length(cfg.subjects)
    cfg.rotcode=['B','C'];
    cfg.subject=cfg.subjects{m};
    
    %% Select ROIs per subject
    switch cfg.experiment{1}
        case 'SVM Gijs'
            switch cfg.subject
                % Gijs Unwanted/Wanted Subject:
                case 'LD';          cfg.ROIs={'V1','V2','V3','V3A','V4D','V7','VP','MT','AIPS','PIPS','V4V','FEF','SEF'};
                case 'MJ';          cfg.ROIs={'V1','V2','V3','V3A','V4V','V7','VP','MT','AIPS','PIPS','FEF','SEF'};
                case 'MM';          cfg.ROIs={'V1','V2','V3','V3A','V4V','V7','VP','MT','AIPS','PIPS'};
                case 'TA';          cfg.ROIs={'V1','V2','V3','V3A','V4V','V7','VP','MT','AIPS','PIPS'};
                case 'TJ';          cfg.ROIs={'V1','V2','V3','V3A','V4V','V7','VP','MT','AIPS','PIPS'};
                case 'TK';          cfg.ROIs={'V1','V2','V3','V3A','V4V','V7','VP','MT', 'AIPS', 'PIPS', 'FEF', 'SEF'};
                case 'TJFCD';       cfg.ROIs={'V1','V2','V3','V3A','V4V','V7','VP','MT','AIPS','PIPS'}; cfg.rotcode=['g','h'];
            end
        case 'SVM Bham'
            switch cfg.subject
                % SVM Bham subjects:
                case 'AM';          cfg.ROIs={'V1', 'V2', 'V3', 'VP', 'V3a', 'V4V', 'V7', 'KO', 'LO', 'pFs', 'MT', 'MST', 'MTplus','AdaptA','AdaptP','AdaptL','MEM','ALTV','ALTA','VIPS','AIPS'};
                case 'HB';          cfg.ROIs={'V1', 'V2', 'V3', 'VP', 'V3a', 'V4V', 'V7', 'KO', 'LO', 'pFs', 'MT', 'MST', 'MTplus','AdaptA','AdaptP','MEM','STAB','ALTV','ALTD','ALTA','VIPS','AIPS'};
                case 'JZ';          cfg.ROIs={'V1', 'V2', 'V3', 'VP', 'V3a', 'V4V', 'V7', 'KO', 'LO', 'pFs', 'MT', 'MST', 'MTplus','AdaptA','AdaptP','AdaptL','MEM','STAB','ALTV','ALTA','VIPS','AIPS','extra'};
                case 'MJ';          cfg.ROIs={'V1', 'V2', 'V3', 'VP', 'V3A', 'V4V', 'V7', 'V8', 'KO', 'IPS', 'MTplus','AdaptA','AdaptL','STAB','ALTV','ALTD','ALTA','VIPS','AIPS'};
                case 'MP';          cfg.ROIs={'V1', 'V2', 'V3', 'VP', 'V3A', 'V4V', 'V7', 'KO', 'LO', 'pFs', 'MT', 'MST', 'MTplus','AdaptA','AdaptP','AdaptL','STAB','ALTV','ALTD','ALTA','VIPS','AIPS','extra'};
                case 'MR';          cfg.ROIs={'V1', 'V2', 'V3', 'VP', 'V3A', 'V4V', 'V7', 'V8', 'KO', 'IPS', 'MTplus','AdaptL','MEM','ALTV','ALTA','VIPS','AIPS','extra'};
                case 'RW';          cfg.ROIs={'V1', 'V2', 'V3', 'VP', 'V3A', 'V4V', 'V7', 'V8', 'KO', 'IPS', 'MTplus','AdaptA','AdaptP','AdaptL','MEM','STAB','ALTV','ALTD','ALTA','VIPS','AIPS'};
            end
        case 'Test'
            switch cfg.subject
                % Test
                case 'TB';        cfg.ROIs={'Test'};
                case 'MJ';        cfg.ROIs={'V1', 'V2', 'V3', 'VP', 'V3A', 'V4V', 'V7', 'V8', 'KO', 'IPS', 'MTplus','AdaptA','AdaptL','STAB','ALTV','ALTD','ALTA','VIPS','AIPS'};
                otherwise
                    cfg.ROIs={'V3A','V1','V2','MT','V4V','VP','V7','V3'};
            end
    end
    
    %remove these 2 lines: 
    %disp ('only calculating a subset of all ROIs');
    %cfg.ROIs = {'V3A','PIPS','AIPS','MT','V2','V1'};
    %cfg.ROIs = {'VP'};
    
    %% Create SWM switches
    disp(' ');
    disp('Creating timecourses from log files');
    [data,cfg]=CreatetimecoursesBatch(cfg);
    
    %% Loop SVM over Subjects
    %  Removing data from Shell input will load the Timecourse files
    switch cfg.shuffletimecourses{1}
        case 'yes';  cfg.shuffleplaceholder = [true];
        case 'no';   cfg.shuffleplaceholder = [false];
        case 'both'; cfg.shuffleplaceholder = [true, false];
    end
    for n = cfg.shuffleplaceholder;
        cfg.shuffletimecourse = n;
        if cfg.shuffletimecourse == 1;    savename = ['ROIs ' datestr(clock,1) ' shuffled'];
        else                              savename = ['ROIs ' datestr(clock,1)]; end;
        clear results
        
        switch cfg.parallelcomputing
            case 0
                for k=1:length(cfg.ROIs)
                    start=tic;
                    results{k}=SVMShell(cfg.ROIs{k},data,cfg,k);
                    disp(['Saving to file: ' cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM Results/'  savename '.mat']);
                    save([cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM Results/' savename '.mat' ])
                    
                    toc(start)
                end;
                
            case 1
                matlabpool;
                parfor k=1:length(ROIs)
                    results{k}=SVMShell(ROIs{k},data,cfg,k);
                end
                matlabpool close
        end
        save([cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM Results/' savename '.mat' ]);
        disp(['saving to file: ' cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM Results/'  savename '.mat']);
    end
end
end

function [data,cfg]=CreatetimecoursesBatch(cfg)
%% CreatetimecoursesBatch
switch cfg.analysistype{1}
    case {'Blocked','Blockedfewfeatures'}
        
        cfg.NoBlocks    = 96;
        cfg.Restlength  = 8;
        cfg.Tasklength  = 9;
        cfg.Skipvols    = [cfg.Restlength+1 cfg.Restlength];
        cfg.BlockLength = cfg.Restlength+cfg.Tasklength;
        
        
        data.Memory     =     ([cfg.datadir cfg.experiment{1} '/' cfg.subject '/' 'LogFiles/Memory Log.dat'])';
        data.Perception = load([cfg.datadir cfg.experiment{1} '/' cfg.subject '/' 'LogFiles/Perception Log.dat'])';
        data.Halves     = [-ones(size(data.Memory).*[0.5,1]); ones(size(data.Memory).*[0.5,1])];
        
        data.Memory     = data.Memory(:)';
        data.Perception = data.Perception(:)';
        data.Halves     = data.Halves(:)';
        
        data.Memorysigned(data.Memory<37 & data.Memory>=0)      = -1;
        data.Memorysigned(data.Memory>=37)       = 1; 
        data.Memorysigned(data.Memory<0)         = 0;
        data.Perceptionsigned(data.Perception<=6 & data.Perception>=-6) = 0;
        data.Perceptionsigned(data.Perception>6) = 1;
        data.Perceptionsigned(data.Perception<-6)= -1;
        
        data.Memorysigned2    = repmat(data.Memorysigned,cfg.BlockLength,1);
        data.Perceptionsigned2= repmat(data.Perceptionsigned,cfg.BlockLength,1);
        data.Halves2          = repmat(data.Halves,cfg.BlockLength,1);
        
        SVMplotswitches(2,data,fieldnames(data),cfg)
        
    case {'FMr matrix','Event related','Multiclass','Beta','BetaAVG','BetaState'}
        for runno=1:cfg.NoLogs
            cfg.runno=runno;
            [data2{runno},cfg] = Createtimecourses(cfg);
        end
        
        data=SVMstructjoin(data2);
        disp(' ');
        
        if (cfg.plottimecourses == 1 && strcmp(cfg.analysistype{1},'Multiclass'))
            handle = 3;
            figure(handle);
            for runno=1:cfg.NoLogs
            
                ypos=floor((runno-1)/4);
                xpos=mod((runno-1),4);
                
                subplot('Position',[xpos/4 0.8-ypos/5 0.25 0.20]);
                hold off
                h=pie([data2{runno}.classSize(1) data2{runno}.classSize(2) data2{runno}.classSize(3) data2{runno}.classSize(4) data2{runno}.classSize(5)]);
                set(h(1),'FaceColor',[1 1 1]);
                set(h(3),'FaceColor',cfg.LLcol);
                set(h(5),'FaceColor',cfg.RLcol);
                set(h(7),'FaceColor',cfg.LRcol);
                set(h(9),'FaceColor',cfg.RRcol);
            end

            subplot('Position',[0 0 1 0.4]);
            hold off
            h=pie([sum(data.classSize(1:5:end)) sum(data.classSize(2:5:end)) sum(data.classSize(3:5:end)) sum(data.classSize(4:5:end)) sum(data.classSize(5:5:end))]);
            set(h(1),'FaceColor',[1 1 1]);
            set(h(3),'FaceColor',cfg.LLcol);
            set(h(5),'FaceColor',cfg.RLcol);
            set(h(7),'FaceColor',cfg.LRcol);
            set(h(9),'FaceColor',cfg.RRcol);
            
            [g gg ggg]=mkdir([cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM results/figures ' cfg.timestarted]);
            disp(['saving timecourseplot: ' cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM results/figures ' cfg.timestarted '/timecourse all sessions.png'])
            saveas(handle,[cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM results/figures ' cfg.timestarted '/timecourse all sessions.png'])

            close(handle);
        end
           
    case {'Test'}
end
end

function data=SVMstructjoin(structs)
%% SVMStructjoin

fns=fieldnames(structs{1});

for m=1:length(fns)
    for n=1:length(structs)
        dim{m,n}=size(structs{n}.(char(fns(m))));
    end
    %data.(char(fns(m)))=zeros(dim{m,n}(1),length(structs)*dim{m,n}(2));
end

for m=1:length(fns)
    for n=1:length(structs)
        data.(char(fns(m)))(1:dim{m,n}(1),((dim{m,n}(2)*(n-1)+1):(dim{m,n}(2)*n)))=...
            structs{n}.(char(fns(m)));
    end
end

end

function [data,cfg]=Createtimecourses(cfg)
%% Createtimecourses
cfg.FlipLength = 200/cfg.SampleRate;

% Read log
[statelist,tasklist,cfg]    = SVMlogread(cfg);                       % Reads Logs in subject/logfiles map
cfg.PrefixLength            = (cfg.onset*1000)/cfg.SampleRate;  % Time it takes to reach steady state (e.g. no onset HRF)
data.TimeCourseLength       = round(cfg.ExpEnd/cfg.SampleRate); %
tasklist(1,2)               = tasklist(1,2)+1;                  % Matlab lists start at 1

% Create Wanted Timecourse (=task)

data.WantedTimecourse       = zeros(data.TimeCourseLength,1);
data.WantedTimecourse(:)    = -tasklist(1,1); % fill in first task (in case task at t=0 is not given

for i  = 1:max(size(tasklist))
    task = tasklist(i,1);
    taskstart = tasklist(i,2);
    if taskstart==0; taskstart=1; end;
    data.WantedTimecourse(taskstart:end) = task;
end

data.WantedTimecourse(1:cfg.PrefixLength) = data.WantedTimecourse(cfg.PrefixLength+1);

% Create State Timecourse

data.StateTimecourse = zeros(data.TimeCourseLength,1);
data.StateTimecourse(:) =  -statelist(1,1); %fill in first percept
for i=1:max(size(statelist))
    state = statelist(i,1);
    statestart = statelist(i,2);
    if statestart==0; statestart=1; end;
    data.StateTimecourse(statestart:end) = state;
end

data.StateTimecourse(1:cfg.PrefixLength) = data.StateTimecourse(cfg.PrefixLength+1);

% Create HoldSwitch Timecourse

data.HoldSwitchTimecourse               = zeros(data.TimeCourseLength,1);
data.HoldSwitchTimecourse(data.StateTimecourse==data.WantedTimecourse) = cfg.Holdcode;
data.HoldSwitchTimecourse(data.StateTimecourse~=data.WantedTimecourse) = -cfg.Holdcode;

    
switch cfg.analysistype{1}
    
    case {'FMr matrix','BetaAVG','BetaState'}
        cfg.originaltimecourses = {'StateTimecourse','WantedTimecourse','HoldSwitchTimecourse','FlipTimecourse','WantedflipTimecourse','UnwantedflipTimecourse','FlipToCCWTimecourse','FlipToCWTimecourse'};
          
        
        %% Create Flip timecourses (old version)
        [rows,cols] = size(data.StateTimecourse);
        for n=4:8; data.(cfg.originaltimecourses{n}) = zeros(rows,cols); end
        
        HalfFlipLength1 = floor(cfg.FlipLength/2);
        HalfFlipLength2 = cfg.FlipLength-HalfFlipLength1;
        
        for i=1:rows-1
            if data.StateTimecourse(i) ~= data.StateTimecourse(i+1)
                if i <= HalfFlipLength1
                    a = HalfFlipLength1+1-i;%correction to avoid references to negative timepoints
                else
                    a=0;
                end
                data.FlipTimecourse((i-HalfFlipLength1+a):(i+HalfFlipLength2-1+a)) = 1;
                if data.WantedTimecourse(i) == data.StateTimecourse(i)
                    data.UnwantedflipTimecourse((i-HalfFlipLength1):(i+HalfFlipLength2-1)) = 1;
                else
                    data.WantedflipTimecourse((i-HalfFlipLength1):(i+HalfFlipLength2-1)) = 1;
                end
                if data.StateTimecourse(i+1) == 1
                    data.FlipToCWTimecourse((i-HalfFlipLength1):(i+HalfFlipLength2-1)) = 1;
                else
                    data.FlipToCCWTimecourse((i-HalfFlipLength1):(i+HalfFlipLength2-1)) = 1;
                end
            end
        end
        
        SVMplotswitches(2,data,cfg.originaltimecourses,cfg);
        
        for n=1:length(cfg.originaltimecourses)
            [data.(['hrfbin' cfg.originaltimecourses{n}]),data.(['hrf' cfg.originaltimecourses{n}])] = SVMhrf(data.(cfg.originaltimecourses{n}),cfg);
        end
        
        for n=1:length(cfg.originaltimecourses); toplot{n}=['hrf' cfg.originaltimecourses{n}]; end
        SVMplotswitches(3,data,toplot,cfg);
        
        for n=1:length(cfg.originaltimecourses); toplot{n}=['hrfbin' cfg.originaltimecourses{n}]; end
        SVMplotswitches(4,data,toplot,cfg);
     
    case 'Multiclass'
        %data = createTrainingQualification(data,cfg);
        data = createTrainingQualificationMulticlass(data,cfg);
        
        cfg.originaltimecourses = {'StateTimecourse','WantedTimecourse','HoldSwitchTimecourse'};
        Stepsize = round((length(data.(cfg.originaltimecourses{1})))/cfg.NoVol); % convert x-axis to cfg.NoVol
        for i=1:size(cfg.originaltimecourses,2)
            data.(['org' cfg.originaltimecourses{i}])=data.(cfg.originaltimecourses{i})(Stepsize/2:Stepsize:end);
        end
    case {'Event related','Beta'}
        %% Create FlipEventRelated timecourses
        
        for n=1:length(cfg.timecourses); data.(cfg.timecourses{n}) = -ones(1,cfg.NoVol); end
        
        Stepsize = round((data.TimeCourseLength)/cfg.NoVol);%convert x-axis to cfg.NoVol
        
        data.FlipERUnwantedTimecourse = zeros(cfg.NoVol,1);
        data.FlipERWantedTimecourse = zeros(cfg.NoVol,1);
        data.FlipERToCWTimecourse = zeros(cfg.NoVol,1);
        data.FlipERToCCWTimecourse = zeros(cfg.NoVol,1);
        data.FlipERWillToCWTimecourse = zeros(cfg.NoVol,1);
        data.FlipERWillToCCWTimecourse = zeros(cfg.NoVol,1);
        data.FlipERTimecourse = zeros(cfg.NoVol,1);
        data.FlipERWillTimecourse = zeros(cfg.NoVol,1);
        
        for i=1:data.TimeCourseLength-1
            if data.StateTimecourse(i) ~= data.StateTimecourse(i+1)
                fliptime = ceil((i+1)/Stepsize); %determine in which volume the flip took place
                if fliptime <= cfg.NoVol
                    data.FlipERTimecourse(fliptime) = 1;
                    if data.WantedTimecourse(i) == data.StateTimecourse(i)
                        data.FlipERUnwantedTimecourse(fliptime) = 1;
                    else
                        data.FlipERWantedTimecourse(fliptime) = 1;
                    end
                    if data.StateTimecourse(i+1) == 1
                        data.FlipERToCWTimecourse(fliptime) = 1;
                    else
                        data.FlipERToCCWTimecourse(fliptime) = 1;
                    end
                end
            end
            
            if data.WantedTimecourse(i) ~= data.WantedTimecourse(i+1)
                fliptime = ceil((i+1)/Stepsize);%determine in which volume the flip took place
                if fliptime <= cfg.NoVol
                    data.FlipERWillTimecourse(fliptime) = 1;
                    if data.WantedTimecourse(i+1) == 1
                        data.FlipERWillToCWTimecourse(fliptime) = 1;
                    else
                        data.FlipERWillToCCWTimecourse(fliptime) = 1;
                    end
                end
            end
        end
        
        for n=1:length(cfg.timecourses)
            data.(cfg.timecourses{n}) = (data.(cfg.timecourses{n})(1:cfg.NoVol));
        end
        
        data.WantedvsUnwantedTimecourse = data.FlipERWantedTimecourse - data.FlipERUnwantedTimecourse;
        data.WillCWvsCCWTimecourse = data.FlipERWillToCWTimecourse - data.FlipERWillToCCWTimecourse;
        
        SVMplotswitches(4,data,cfg.timecourses,cfg)
        
end
end

function data = createTrainingQualificationMulticlass(data,cfg)
        %L=CW, R=CCW. LR->wanted left, state right 
        tmp.LL = data.WantedTimecourse==-cfg.CCWcode & data.StateTimecourse==-cfg.CCWcode;
        tmp.RL = data.WantedTimecourse== cfg.CCWcode & data.StateTimecourse==-cfg.CCWcode;
        tmp.LR = data.WantedTimecourse==-cfg.CCWcode & data.StateTimecourse== cfg.CCWcode;
        tmp.RR = data.WantedTimecourse== cfg.CCWcode & data.StateTimecourse== cfg.CCWcode;

        p = [6 16 1 1 6 0 32];      % Is default
        hrf = spm_hrf_backup(1/10, p);

        hrfLLcourse=conv(double(tmp.LL),hrf);
        hrfLLcourse=hrfLLcourse(1:length(tmp.LL));
        hrfRLcourse=conv(double(tmp.RL),hrf);
        hrfRLcourse=hrfRLcourse(1:length(tmp.RL));
        hrfLRcourse=conv(double(tmp.LR),hrf);
        hrfLRcourse=hrfLRcourse(1:length(tmp.LR));
        hrfRRcourse=conv(double(tmp.RR),hrf);
        hrfRRcourse=hrfRRcourse(1:length(tmp.RR));

        tmp.Stepsize = round((length(tmp.LL))/cfg.NoVol); % convert x-axis to cfg.NoVol
        tmp.rescalehrfLLcourse=hrfLLcourse(tmp.Stepsize/2:tmp.Stepsize:end);
        tmp.rescalehrfRLcourse=hrfRLcourse(tmp.Stepsize/2:tmp.Stepsize:end);
        tmp.rescalehrfLRcourse=hrfLRcourse(tmp.Stepsize/2:tmp.Stepsize:end);
        tmp.rescalehrfRRcourse=hrfRRcourse(tmp.Stepsize/2:tmp.Stepsize:end);

        %select only the training points that are significantly in one
        %class. i.e. have a bold signal that is higher than the other
        %signals with a certain threshold.
        data.LLselection = (tmp.rescalehrfLLcourse - max([tmp.rescalehrfRLcourse tmp.rescalehrfLRcourse tmp.rescalehrfRRcourse],[],2)) > 0.2;
        data.RLselection = (tmp.rescalehrfRLcourse - max([tmp.rescalehrfLLcourse tmp.rescalehrfLRcourse tmp.rescalehrfRRcourse],[],2)) > 0.2;
        data.LRselection = (tmp.rescalehrfLRcourse - max([tmp.rescalehrfRLcourse tmp.rescalehrfLLcourse tmp.rescalehrfRRcourse],[],2)) > 0.2;
        data.RRselection = (tmp.rescalehrfRRcourse - max([tmp.rescalehrfRLcourse tmp.rescalehrfLRcourse tmp.rescalehrfLLcourse],[],2)) > 0.2;
            
        %disp('data.LLselection = ');
        %disp(data.LLselection);
        
        data.WantedTrainingQualification = (data.LLselection + data.LRselection) - (data.RLselection + data.RRselection) * -cfg.CCWcode;
        data.StateTrainingQualification =  (data.LLselection + data.RLselection) - (data.LRselection + data.RRselection) * -cfg.CCWcode;
                
        data.LLvsLRTimecourse = data.LLselection - data.LRselection;
        data.RLvsRRTimecourse = data.RLselection - data.RRselection;
        data.LLvsRLTimecourse = data.LLselection - data.RLselection;
        data.LRvsRRTimecourse = data.LRselection - data.RRselection;
        data.LLvsRRTimecourse = data.LLselection - data.RRselection;
        data.LRvsRLTimecourse = data.LRselection - data.RLselection;
        
        %data.classSize(1)=(sum((data.StateTrainingQualification .* data.WantedTrainingQualification)==0));
        %data.classSize(2)=(sum((data.StateTrainingQualification==-cfg.CCWcode) .* (data.WantedTrainingQualification==-cfg.CCWcode)));
        %data.classSize(3)=(sum((data.StateTrainingQualification==-cfg.CCWcode) .* (data.WantedTrainingQualification== cfg.CCWcode)));
        %data.classSize(4)=(sum((data.StateTrainingQualification== cfg.CCWcode) .* (data.WantedTrainingQualification==-cfg.CCWcode)));
        %data.classSize(5)=(sum((data.StateTrainingQualification== cfg.CCWcode) .* (data.WantedTrainingQualification== cfg.CCWcode)));
        data.classSize(1) = sum(~data.LLselection & ~data.LRselection & ~data.RLselection & ~data.RRselection);
        data.classSize(2) = sum(data.LLselection);
        data.classSize(3) = sum(data.LRselection);
        data.classSize(4) = sum(data.RLselection);
        data.classSize(5) = sum(data.RRselection);
        
        disp (['class 0   : ' int2str(data.classSize(1))])
        disp (['class 1 LL: ' int2str(data.classSize(2))])
        disp (['class 2 RL: ' int2str(data.classSize(3))])
        disp (['class 3 LR: ' int2str(data.classSize(4))])
        disp (['class 4 RR: ' int2str(data.classSize(5))])
           
        if (cfg.plottimecourses == 1)
            plotTrainingQualificationMulticlass(2,data,tmp,cfg);
        end
end

function plotTrainingQualificationMulticlass(handle,data,tmp,cfg)
    figure(handle);
    
    set(gcf,'DefaultAxesColorOrder',[cfg.LLcol;cfg.RLcol;cfg.LRcol;cfg.RRcol;0 0 1;1 1 0]);
    
    subplot('Position',[0.05 0.6 0.7 0.35]);
    plot([16+tmp.LL(tmp.Stepsize/2:tmp.Stepsize:end) 13+tmp.RL(tmp.Stepsize/2:tmp.Stepsize:end) 10+tmp.LR(tmp.Stepsize/2:tmp.Stepsize:end) 7+tmp.RR(tmp.Stepsize/2:tmp.Stepsize:end) 4+data.WantedTimecourse(tmp.Stepsize/2:tmp.Stepsize:end) 1+data.StateTimecourse(tmp.Stepsize/2:tmp.Stepsize:end)]);   
    legend('LL','RL','LR','RR','Wanted','State','Location','NorthOutside','Orientation','horizontal');
    
    subplot('Position',[0.05 0.3 0.7 0.25]);
    hold off
    plot([tmp.rescalehrfLLcourse tmp.rescalehrfRLcourse tmp.rescalehrfLRcourse tmp.rescalehrfRRcourse]);

    subplot('Position',[0.05 0.05 0.7 0.15]);
    hold off
    colormap([cfg.RRcol; 1 1 1; cfg.RLcol; 1 1 1; cfg.LRcol; 1 1 1; cfg.LLcol]);
    image(4 + (data.StateTrainingQualification + 2*data.WantedTrainingQualification)');
    
    subplot('Position',[0.75 0.05 .25 0.25]);
    hold off
    h=pie([data.classSize(1) data.classSize(2) data.classSize(3) data.classSize(4) data.classSize(5)]);
    set(h(1),'FaceColor',[1 1 1]);
    set(h(3),'FaceColor',cfg.LLcol);
    set(h(5),'FaceColor',cfg.RLcol);
    set(h(7),'FaceColor',cfg.LRcol);
    set(h(9),'FaceColor',cfg.RRcol);

    [g gg ggg]=mkdir([cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM results/figures ' cfg.timestarted]);
    disp(['saving timecourseplot: ' cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM results/figures ' cfg.timestarted '/timecourse session ' num2str(cfg.runno) '.png'])
    saveas(handle,[cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM results/figures ' cfg.timestarted '/timecourse session ' num2str(cfg.runno) '.png'])
    
    close(handle);
end

function data = createTrainingQualification(data,cfg)
        [data.hrfSeparateStateTimecourse1 data.hrfSeparateStateTimecourse2] = SVMSeparatehrf(data.StateTimecourse,cfg);
        [data.hrfSeparateWantedTimecourse1 data.hrfSeparateWantedTimecourse2] = SVMSeparatehrf(data.WantedTimecourse,cfg);
        
        %todo: 0.4 in cfg variabele zetten
        %timecourse1=CCW, timecourse2=CW, dus l=cw, r=ccw;
        sl = (data.hrfSeparateStateTimecourse2 - data.hrfSeparateStateTimecourse1) > 0.1;
        sr = (data.hrfSeparateStateTimecourse1 - data.hrfSeparateStateTimecourse2) > 0.1;
        wl = (data.hrfSeparateWantedTimecourse2 - data.hrfSeparateWantedTimecourse1) > 0.1;
        wr = (data.hrfSeparateWantedTimecourse1 - data.hrfSeparateWantedTimecourse2) > 0.1;
        
        data.StateTrainingQualification = sl-sr;
        data.WantedTrainingQualification = wl-wr;
        
        %deze kan je gebruiken voor de one vs all.
        %data.TrainWantedCCWStateCCW = -1 + 2 * ((data.StateTrainingQualification== cfg.CCWcode) .* (data.WantedTrainingQualification== cfg.CCWcode));
        %data.TrainWantedCCWStateCW  = -1 + 2 * ((data.StateTrainingQualification==-cfg.CCWcode) .* (data.WantedTrainingQualification== cfg.CCWcode));
        %data.TrainWantedCWStateCCW  = -1 + 2 * ((data.StateTrainingQualification== cfg.CCWcode) .* (data.WantedTrainingQualification==-cfg.CCWcode));
        %data.TrainWantedCWStateCW   = -1 + 2 * ((data.StateTrainingQualification==-cfg.CCWcode) .* (data.WantedTrainingQualification==-cfg.CCWcode));
        
        %deze voor de all vs all.
        %data.TrainWantedCCWStateCCW = wr .* sr; %RR %((data.StateTrainingQualification== cfg.CCWcode) .* (data.WantedTrainingQualification== cfg.CCWcode));
        %data.TrainWantedCCWStateCW  = wr .* sl; %RL %((data.StateTrainingQualification==-cfg.CCWcode) .* (data.WantedTrainingQualification== cfg.CCWcode));
        %data.TrainWantedCWStateCCW  = wl .* sr; %LR %((data.StateTrainingQualification== cfg.CCWcode) .* (data.WantedTrainingQualification==-cfg.CCWcode));
        %data.TrainWantedCWStateCW   = wl .* sl; %LL %((data.StateTrainingQualification==-cfg.CCWcode) .* (data.WantedTrainingQualification==-cfg.CCWcode));
        
        data.LLvsLRTimecourse = (wl & sl) - (wl & sr);
        data.RLvsRRTimecourse = (wr & sl) - (wr & sr);
        data.LLvsRLTimecourse = (wl & sl) - (wr & sl);
        data.LRvsRRTimecourse = (wl & sr) - (wr & sr);
        data.LLvsRRTimecourse = (wl & sl) - (wr & sr);
        data.LRvsRLTimecourse = (wl & sr) - (wr & sl);
        
        disp (['0: ' int2str((sum((data.StateTrainingQualification .* data.WantedTrainingQualification)==0)))])
        disp (['1: ' int2str((sum((data.StateTrainingQualification==-cfg.CCWcode) .* (data.WantedTrainingQualification==-cfg.CCWcode))))])
        disp (['2: ' int2str((sum((data.StateTrainingQualification==-cfg.CCWcode) .* (data.WantedTrainingQualification== cfg.CCWcode))))])
        disp (['3: ' int2str((sum((data.StateTrainingQualification== cfg.CCWcode) .* (data.WantedTrainingQualification==-cfg.CCWcode))))])
        disp (['4: ' int2str((sum((data.StateTrainingQualification== cfg.CCWcode) .* (data.WantedTrainingQualification== cfg.CCWcode))))])
               
        %werkt niet meer omdat 1 timecourse al regescaled is en de ander
        %niet..
        %plotTrainingQualification(2,data,cfg);
end

function plotTrainingQualification(handle,data,cfg)
figure(handle);

colormap([0.6 0.6 1; 1 1 1; 0.6 1 0.6]);

subplot(3,1,1);
hold on
image(2 + data.StateTrainingQualification');
plot([data.hrfSeparateStateTimecourse1 data.hrfSeparateStateTimecourse2 data.StateTimecourse])
hold off

subplot(3,1,2);
hold on
image(2 + data.WantedTrainingQualification');
plot([data.hrfSeparateWantedTimecourse1 data.hrfSeparateWantedTimecourse2 data.WantedTimecourse])
hold off

subplot(3,1,3);
image(2 + [data.StateTrainingQualification data.WantedTrainingQualification]');

end

function [statelist,tasklist,cfg]=SVMlogread(cfg)
%% SVMlogread

fname = [cfg.datadir cfg.experiment{1} '/' cfg.subject '/LogFiles/' cfg.subject '_BistableSVMIntention-logfile' num2str(cfg.runno) '.dat'];

disp(fname);
fid = fopen(fname,'r');

for i = 1:cfg.NoOfEmptyLines+1
    tline = fgetl(fid); %empty line
end

cfg.TrialDur=[];
while ~isempty(sscanf(tline(1),'%d'))
    cfg.TrialDur(end+1,:) = sscanf(tline, '%d %d');
    tline = fgetl(fid);
end

curtime = 0;
ResponseNo = 0;
TrialNo = 0;
TimeCheck = [];
while 1
    if ~ischar(tline), break, end
    a1 = strfind(tline,['Received : ' cfg.rotcode(1)]);
    a2 = strfind(tline,['Received : ' cfg.rotcode(2)]);
    a3 = strfind(tline,'See Right at time');
    a4 = strfind(tline,'See Left at time');
    a5 = strfind(tline,'Experiment Lasted');
    if (~isempty(a1))
        ResponseNo = ResponseNo + 1;
        [f1,f2] = strtok(tline,':'); %splits the string in 2 at indicated position
        [f3,f4] = strtok(f2(2:end),':');
        curtime = curtime + str2num(f4(2:end));
        triggertimeL = curtime - cfg.RespDelayCor;
        statelist(ResponseNo,1) = -cfg.CCWcode;
        statelist(ResponseNo,2) = triggertimeL;
        TimeCheck = [TimeCheck curtime];
    end
    if (~isempty(a2))
        ResponseNo = ResponseNo + 1;
        [f1,f2] = strtok(tline,':');
        [f3,f4] = strtok(f2(2:end),':');
        curtime = curtime + str2num(f4(2:end));
        triggertimeR = curtime - cfg.RespDelayCor;
        statelist(ResponseNo,1) = cfg.CCWcode;
        statelist(ResponseNo,2) = triggertimeR;
        TimeCheck = [TimeCheck curtime];
    end
    if (~isempty(a3))
        TrialNo = TrialNo + 1;
        [f1,f2] = strtok(tline,':');
        TrialStart = str2num(f2(2:end));
        tasklist(TrialNo,1) = cfg.CCWcode;
        tasklist(TrialNo,2) = TrialStart;
        curtime = TrialStart;
        TimeCheck = [TimeCheck curtime];
    end
    if (~isempty(a4))
        TrialNo = TrialNo + 1;
        [f1,f2] = strtok(tline,':');
        TrialStart = str2num(f2(2:end));
        tasklist(TrialNo,1) = -cfg.CCWcode;
        tasklist(TrialNo,2) = TrialStart;
        curtime = TrialStart;
        TimeCheck = [TimeCheck curtime];
    end
    if (~isempty(a5))
        [f1,f2] = strtok(tline,':');
        cfg.ExpEnd = str2num(f2(2:end));
        curtime = cfg.ExpEnd;
        TimeCheck = [TimeCheck curtime];
    end
    tline = fgetl(fid);
end

if statelist(2,end)==0
    disp('Change cfg.rotcode to suit the experiment!')
end;

% Lower SampleRate
%% double check the timeshift thing!
statelist(:,2) = round((statelist(:,2) + cfg.TimeShift)/cfg.SampleRate);
tasklist(:,2) = round((tasklist(:,2) + cfg.TimeShift)/cfg.SampleRate);
tasklist((tasklist(:,2) <= 0), 2) = 0;          %remove times below zero
statelist((statelist(:,2) <= 0), 2) = 0;
cfg.NoVol=round(sum(cfg.TrialDur(:,1))/cfg.Tr);

fclose(fid);
end

function SVMplotswitches(handle,data,timecourses,cfg)
%% Plot Switches
if (cfg.plottimecourses == 1)
    pos=mod(handle,8);
    figure(handle);
    if ~isequal(get(handle,'WindowStyle'),'docked')
        set(handle,'Position',[pos*100 0 400 700]);
    end
    set(handle, 'Name', ['CCW = ' num2str(cfg.CCWcode) '   Hold = ' num2str(cfg.Holdcode)]);
    hold on
    for n=1:length(timecourses)
        hold off
        Handle = subplot(length(timecourses),1,n);
        imagesc(data.(timecourses{n}));
        colormap(gray);
        %         whitebg('k');
        
        set(handle, 'Color', [0.75 0.75 0.75]);
        
        box off;
        title(timecourses{n});
        %hold off
    end;
end

end

function [binhrfcourse,hrfcourse] = SVMhrf(timecourse,cfg)
%% hrfcourse
if nargin<3; verbose=0; else verbose=1; end
% Plots only when handle is filled in.
correction=0.0000000001;    % Was in the original code, probably enables signing if =0;
p = [6 16 1 1 6 0 32];      % Is default
hrf = spm_hrf_backup(1/10, p);
hrfcourse=conv(double(timecourse==max(timecourse)),hrf)-conv(double(timecourse==min(timecourse)),hrf);
hrfcourse=hrfcourse(1:length(timecourse));

Stepsize = round((length(timecourse))/cfg.NoVol); % convert x-axis to cfg.NoVol
binhrfcourse=sign(hrfcourse(Stepsize/2:Stepsize:end)-correction);

end

function [rescalehrfcourse1,rescalehrfcourse2] = SVMSeparatehrf(timecourse,cfg)
%% hrfcourse
p = [6 16 1 1 6 0 32];      % Is default
hrf = spm_hrf_backup(1/10, p);

hrfcourse1=conv(double(timecourse==cfg.CCWcode),hrf);
hrfcourse1=hrfcourse1(1:length(timecourse));
hrfcourse2=conv(double(timecourse==-cfg.CCWcode),hrf);
hrfcourse2=hrfcourse2(1:length(timecourse));

Stepsize = round((length(timecourse))/cfg.NoVol); % convert x-axis to cfg.NoVol
rescalehrfcourse1=hrfcourse1(Stepsize/2:Stepsize:end);
rescalehrfcourse2=hrfcourse2(Stepsize/2:Stepsize:end);
end

function results=SVMShell(ROI,data,cfg,id,varargin)
%% SVMShell - Import data and make useable
% Everything stored in results will be saved to file each run in  a cell
% results{ROInumber};

numoptargs = nargin -4;
optargs = varargin;
results.subject =cfg.subject;

if  numoptargs >= 1
    if rem(numoptargs,2)~=0
        error('Bioinfo:svmtrain:IncorrectNumberOfArguments',...
            'Incorrect number of arguments to %s.',mfilename);
    end
    for n=1:2:numoptargs
        cfg.(optargs{n})=optargs{n+1};
    end
end

loaded=0;
while loaded==0;
    try
        data.id  = num2str(id);
        data.ROI    = ROI;
        name1 = [cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM Data/' 'timecourses_' cfg.subject '-NOSMOOTHING.mat'];
        load(name1,'matrix');
        disp(['loading: ' name1]);
        loaded=1;
    catch ME
        disp(ME.message);
        loaded=0;
    end
end

%%%%%%%%%%%%%%%%% 
%[matrix,data] = generateFakeData(matrix,data,cfg);
%%%%%%%%%%%%%%%%%

%% Create a matrix indicating which volume belongs to which run:
cfg.NoVox         = size(matrix.(data.ROI),1);
cfg.NoVol         = size(matrix.(data.ROI),2);

disp(['Contains ' num2str(cfg.NoVox) ' Voxels, and ' num2str(cfg.NoVol) ' Volumes']);

% Preprocessing
matrix.currentROI = zeros(size(matrix.(data.ROI)));

switch cfg.normalization{1} % Normalization, based on 2008 paper
    case 'mean block, std scan'
        for z=1:size(matrix.(data.ROI),1)
            for y=1:size(matrix.(data.ROI),3)
                if std(matrix.(data.ROI)(z,:))>0
                    matrix.currentROI(z,:,y)= (matrix.(data.ROI)(z,:,y)-mean(matrix.(data.ROI)(z,:,y)))/std(matrix.(data.ROI)(z,:)); else
                    matrix.currentROI(z,:,y)= (matrix.(data.ROI)(z,:,y)-mean(matrix.(data.ROI)(z,:,y)));
                end
            end
        end
        data.maxstd      = 3;
        matrix.currentROI(matrix.currentROI>data.maxstd) = data.maxstd;
        matrix.currentROI(matrix.currentROI<-data.maxstd)= -data.maxstd;
        disp(['Number of outliers found during normalization: ' num2str(sum(abs(matrix.currentROI(:))==data.maxstd)) ' of ' num2str(length(matrix.currentROI(:)))  ' voxels ; Standard Deviation > ' num2str(data.maxstd)])
    case 'mean scan, std scan'
        for z=1:size(matrix.(data.ROI),1)
            for y=1:size(matrix.(data.ROI),3)
                matrix.currentROI(z,:,y)= (matrix.(data.ROI)(z,:,y)-mean(matrix.(data.ROI)(z,:)))/std(matrix.(data.ROI)(z,:));
            end
        end
        data.maxstd      = 3;
        matrix.currentROI(matrix.currentROI>data.maxstd) = data.maxstd;
        matrix.currentROI(matrix.currentROI<-data.maxstd)= -data.maxstd;
        disp(['Number of outliers found during normalization: ' num2str(sum(abs(matrix.currentROI(:))==data.maxstd)) ' of ' num2str(length(matrix.currentROI(:)))  ' voxels ; Standard Deviation > ' num2str(data.maxstd)])
    case 'mean'
        matrix.currentROI= matrix.currentROI - repmat(mean(mean(matrix.currentROI,3),2),[1,size(matrix.currentROI,2),size(matrix.currentROI,3)]);
end

matrix.BandROI   = bandpass(cfg,matrix.currentROI);
matrix.SmoothROI = supersmooth(cfg,matrix.BandROI);

% For testing purposes, creates matrix
% if isequal(cfg.experiment{1},'Test')
% %    matrix.SmoothROI = reshape(repmat(1:size(matrix.SmoothROI,2)*size(matrix.SmoothROI,3)',[size(matrix.SmoothROI,1),1,1]),size(matrix.SmoothROI));
%
%     matrix.SmoothROI = ones(size(matrix.SmoothROI));
%     matrix.SmoothROI(:,:,[2 4]) = 0;
% end


%% Create matrix from events
for p = 1:length(cfg.timecourses);
    
    % Create matrices
    disp(['Creating voxel timecourse matrix for ROI #' data.id ': ' data.ROI])
    switch cfg.analysistype{1}
        case 'Blocked'
            %birmingham, project van maartje met 20s uit, 20 sec
            %intermittent data.
            
            cfg.NoTimepoints  = size(data.Halves,1);
            data.matrix       = zeros(cfg.NoVox,cfg.NoVol,cfg.NoLogs);
            data.relvol       = ones(cfg.NoTimepoints,cfg.NoLogs);  % Determines relavant volumes
            % Correct for missing volumes / high intensity volumes
            for n=1:cfg.NoLogs
                data.matrix(:,1+cfg.fmrvolumestoskip:end,n) = matrix.SmoothROI(:,1+cfg.fmrvolumestoskip:cfg.NoVol,n);
            end;
            
            data.matrix     = data.matrix(:,cfg.Skipvols(1):(end-cfg.Skipvols(2)),:);
            data.matrix     = reshape(data.matrix,[size(data.matrix,1),cfg.BlockLength,cfg.NoBlocks]);
            
            data.actualuseTrain   = data.(cfg.timecourses{p})~=0;
            data.actualuseTest    = data.(cfg.timecourses{p})~=0;
            cfg.NoInstances = cfg.NoBlocks;
            
            for n=1:cfg.NoBlocks
                temp              = data.matrix(:,:,n);
                data.fmrmatrix{n} = temp(:);
            end
        case 'Blockedfewfeatures'
            % zelfde als blocked, maar da alleen de 
            %
            
            tic
            cfg.NoTimepoints= size(data.Halves,1);
            data.matrix     = zeros(cfg.NoVox,cfg.NoVol,cfg.NoLogs);
            data.relvol     = data.(cfg.timecourses{p})~=0;  % Determines relavant volumes
            % Correct for missing volumes / high intensity volumes
            for n=1:cfg.NoLogs
                data.matrix(:,1+cfg.fmrvolumestoskip:end,n) = matrix.SmoothROI(:,1+cfg.fmrvolumestoskip:cfg.NoVol,n);
            end;
            
            hrf             = spm_hrf_backup(cfg.Tr/1000); %Only high intensity volumes are relevant
            event           = [ones(1,cfg.Tasklength) zeros(1,cfg.Restlength)];
            data.hrfcourse  = conv(  double(event),hrf);
            data.hrfcourse  = data.hrfcourse(1:cfg.BlockLength)>0.20*max(data.hrfcourse);
            
            data.matrix     = data.matrix(:,cfg.Skipvols(1):(end-cfg.Skipvols(2)),:);
            data.matrix     = reshape(data.matrix,[size(data.matrix,1),cfg.BlockLength,cfg.NoBlocks]);
            
            data.actualuseTrain  = repmat(data.hrfcourse,1,cfg.NoBlocks) .* data.relvol;
            data.actualuseTest  =  repmat(data.hrfcourse,1,cfg.NoBlocks) .* data.relvol;
            cfg.NoInstances = cfg.NoBlocks;
            
            
            for n=1:cfg.NoBlocks;
                data.fmrmatrix{n} = data.matrix(:,:,n);
            end
            
            disp(['time taken to create matrix = ' num2str(toc) ' seconds']);
            
        case 'FMr matrix'
            
            cfg.NoTimepoints  = size(data.(cfg.timecourses{1}),1);
            data.matrix       = zeros(cfg.NoVox,cfg.NoTimepoints,cfg.NoLogs);
            data.relvol       = zeros(cfg.NoTimepoints,cfg.NoLogs);             % Determines relavant volumes
            % Correct for missing volumes / high intensity volumes
            for n=1:cfg.NoLogs
                data.matrix(:,1+cfg.NoTimepoints-cfg.NoVol+cfg.fmrvolumestoskip:end,n) = matrix.SmoothROI(:,1+cfg.fmrvolumestoskip+(cfg.NoVol-cfg.NoVol):cfg.NoVol,n);
                % Determine if volumes are relavant
                data.relvol(find(data.(cfg.timecourses{1})(:,n)~=data.(cfg.timecourses{1})(1,n),1):end,n)=1;
            end;
            
            cfg.NoInstances           = cfg.NoLogs;
            for n=1:cfg.NoLogs
                data.fmrmatrix{n}          = data.matrix(:,:,n);
            end
            
            data.actualuseTrain              = data.relvol;
            data.actualuseTest               = data.relvol;
        case 'Multiclass'
            
            cfg.NoTimepoints  = size(data.(cfg.timecourses{1}),1);
            data.matrix       = zeros(cfg.NoVox,cfg.NoTimepoints,cfg.NoLogs);
            data.relvol       = zeros(cfg.NoTimepoints,cfg.NoLogs);             % Determines relavant volumes
            % Correct for missing volumes / high intensity volumes
            for n=1:cfg.NoLogs
                data.matrix(:,1+cfg.NoTimepoints-cfg.NoVol+cfg.fmrvolumestoskip:end,n) = matrix.SmoothROI(:,1+cfg.fmrvolumestoskip+(cfg.NoVol-cfg.NoVol):cfg.NoVol,n);
                % Determine if volumes are relavant
                data.relvol(find(data.(cfg.timecourses{1})(:,n)~=data.(cfg.timecourses{1})(1,n),1):end,n)=1;
                data.fmrmatrix{n}          = data.matrix(:,:,n);
            end;
            
            cfg.NoInstances           = cfg.NoLogs;
            
            data.actualuseTrain             = data.relvol .* (data.(cfg.timecourses{p})~=0); %train only on data where the timecourse is +1 or -1, a timecourse value of zero means that it is no significant training data 
            data.actualuseTest              = data.relvol; 
            
        case 'Event related'
            
            cfg.NoTimepoints  = size(data.(cfg.timecourses{1}),1);
            data.matrix       = zeros(cfg.NoVox,cfg.NoTimepoints,cfg.NoLogs);
            data.relvol       = zeros(cfg.NoTimepoints,cfg.NoLogs);             % Determines relavant volumes
            % Correct for missing volumes / high intensity volumes
            for n=1:cfg.NoLogs
                data.matrix(:,1+cfg.NoTimepoints-cfg.NoVol+cfg.fmrvolumestoskip:end,n) = matrix.SmoothROI(:,1+cfg.fmrvolumestoskip+(cfg.NoVol-cfg.NoVol):cfg.NoVol,n);
                % Determine if volumes are relavant
                data.relvol(find(data.(cfg.timecourses{1})(:,n)~=data.(cfg.timecourses{1})(1,n),1):end,n)=1;
            end;
            cfg.NoInstances            = cfg.NoLogs;
            data.fmrmatrix             = cell(1,cfg.NoLogs);
            data.actualuseTrain              = zeros(cfg.NoTimepoints,cfg.NoLogs);
            for n=1:cfg.NoLogs
                
                cfg.runlength(n)           = length(find(data.relvol(:,n)))-cfg.CombineEvents(1)-cfg.CombineEvents(2);
                data.fmrmatrix{n}          = zeros(cfg.NoVox*(cfg.CombineEvents(2)-cfg.CombineEvents(1)+1),cfg.runlength(n));
                data.actualuseTrain(find(data.relvol(:,n),1)+cfg.CombineEvents(1):find(data.relvol(:,n),1,'last')-cfg.CombineEvents(2),n) = 1;
                
                for m=cfg.CombineEvents(1):cfg.CombineEvents(2)
                    data.fmrmatrix{n}(cfg.NoVox*(m-cfg.CombineEvents(1))+1:cfg.NoVox*(m-cfg.CombineEvents(1)+1),:)...
                        = data.matrix(:,cfg.CombineEvents(1)+m:cfg.CombineEvents(2)+m,n);
                end
            end
            data.actualuseTest=data.actualuseTrain;
            
        case {'Beta','BetaAVG','BetaState'}
            
            
            cfg.NoTimepoints  = size(data.(cfg.timecourses{p}),1);
            data.relvol       = ones(2,cfg.NoLogs);             % Determines relavant volumes           
            cfg.NoInstances           = cfg.NoLogs;
            
            data.matrix       = zeros(cfg.NoVox,cfg.NoTimepoints,cfg.NoLogs);

            for n=1:cfg.NoLogs
                data.matrix(:,1+cfg.NoTimepoints-cfg.NoVol+cfg.fmrvolumestoskip:end,n) = ...
                matrix.SmoothROI(:,1+cfg.fmrvolumestoskip+(cfg.NoVol-cfg.NoVol):cfg.NoVol,n);
            end;
            
            data.(['BetaAtypical' cfg.timecourses{p}]) = zeros(2,cfg.NoLogs);
            for n=1:cfg.NoLogs
                if strcmp(cfg.analysistype{1},'BetaAVG')
                    data.fmrmatrix{n}                = findBetasAVG(data.matrix(:,:,n),data.(cfg.timecourses{p})(:,n)',cfg)';
                else
                    data.fmrmatrix{n}                = findBetas2(data.matrix(:,:,n),data.(cfg.timecourses{p})(:,n)',cfg)';    
                end
                data.(['BetaAtypical' cfg.timecourses{p}])(:,n) = [-1 1];
            end
            
            cfg.timecourses{p} = ['BetaAtypical' cfg.timecourses{p}];
            
            if cfg.plottimecourses == 1
                figure(8);
                imagesc(data.fmrmatrix{1});
            end
            
            data.actualuseTrain              = data.relvol;
            data.actualuseTest               = data.relvol;         
    end
    
    
    if cfg.exportmatrix == 1;
        cfg.timecoursestoanalyse=[];
        [xx,yy,zz]=mkdir([cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM Matrices']);
        save(            [cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM Matrices/'  data.ROI ' .mat']);
        disp(['Saving '   cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM Matrices/'  data.ROI ' .mat']);
    end
    %% Create training and test data
    
    disp('Creating training and test data');
    data.runlist=1:cfg.NoInstances;
    
    rand('seed',2);
    
    for n=1:cfg.SVMruns
        cfg.Testruns{n}=data.runlist(randperm(cfg.NoInstances));
        cfg.TestrunsIndex{n}=cfg.Testruns{n}(1:(round(cfg.NoInstances*cfg.trainsetsize)));
    end
    
    if cfg.randomdata == 1
        [d1,d2] = size(cfg.fmrmatrix);
        cfg.fmrmatrix = rand(d1,d2);
        disp('RANDOM DATA!');
    end
    
    for SVMpass = 1:cfg.SVMruns
        
        cfg.Testrunlogical       = zeros(size(data.(cfg.timecourses{p})));
        cfg.Testrunlogical(:,cfg.TestrunsIndex{SVMpass}) = 1;
        
        % Creates actual timecourses using actualuse and testrunlogical
            data.trainmatrix       = SVMjoincell(data.fmrmatrix, data.actualuseTrain.*cfg.Testrunlogical);
            data.testmatrix        = SVMjoincell(data.fmrmatrix, data.actualuseTest.*not(cfg.Testrunlogical));
            data.trainevents       = data.(cfg.timecourses{p})(logical(data.actualuseTrain.*cfg.Testrunlogical));
            data.testevents        = data.(cfg.timecourses{p})(logical(data.actualuseTest.*not(cfg.Testrunlogical)));

          if (cfg.plottimecourses == 1)
            figure(9);
            imagesc(data.trainmatrix);
          end
          
        % Shuffle time courses
        if cfg.shuffletimecourse == 1
            data.trainevents = data.trainevents(randperm(size(data.trainevents,1)),:);
        end
        %% Switch between Machine Learners
        cfg.infotext = [cfg.subject ',' data.ROI ',' char(cfg.timecourses{p}) ',shuffled ' num2str(cfg.shuffletimecourse) ',' cfg.timestarted];
        
        if cfg.shuffletimecourse == 1
            disp(['Running ' cfg.SVMtype{1} ', Shuffled, Machine Learner for ROI thread: ' data.id ', Timecourse: ' (cfg.timecourses{p}) ', ' data.ROI ', Contains: ' num2str(size(data.trainmatrix,1)) ' Features and ' num2str(size(data.testmatrix,2)+size(data.trainmatrix,2)) ' Timepoints']);else
            disp(['Running ' cfg.SVMtype{1} ' Machine Learner for ROI thread: ' data.id ', Timecourse: ' (cfg.timecourses{p}) ', ' data.ROI ', Contains: ' num2str(size(data.trainmatrix,1)) ' Features and ' num2str(size(data.testmatrix,2)+size(data.trainmatrix,2)) ' Timepoints']);
        end
        
        disp('Dimesnsion of data:');
        disp(size(data.fmrmatrix));
        
        switch cfg.SVMtype{1}
            
            case 'winSVM'
                results = SVMDos    (cfg,data,cfg.timecourses{p},SVMpass,results);
                
            case 'Matlab'
                results = SVMMatlab (cfg,data,cfg.timecourses{p},SVMpass,results);
                
            case 'Weka'
                results = SVMDos    (cfg,data,cfg.timecourses{p},SVMpass,results);
                
            case 'Mex'
                error('Not yet implemented');
        end
        
        %% Calculate Accuracies & Plot
        % Syntax -> Results.(variant).data.variable
        % testprediction:   Predicted events during test set
        % normal:           All classes counted
        
        % scaledtestprediction: Same but Number of events = number predicted
        % raw:              Minority class counted only
        % normalized:       Number of events = number predicted (scaled)
        % Chancelevels:     Chance of returning minority class
        results.ROI = ROI;
        
        results.(cfg.timecourses{p}).testevents{SVMpass}  = data.testevents;
        results.(cfg.timecourses{p}).trainevents{SVMpass} = data.trainevents;
        results.minclass                          = 1; % Minority class
        
        % CALCULATE NON-NORMALIZED ACCURACY:
        results.(cfg.timecourses{p}).nodetected(SVMpass)     = sum(sign(results.(cfg.timecourses{p}).predTimecourse{SVMpass}) == results.minclass);
        results.(cfg.timecourses{p}).nopresent(SVMpass)      = sum(data.testevents == results.minclass);
        results.(cfg.timecourses{p}).accuracy.raw(SVMpass)   = sum(results.(cfg.timecourses{p}).testprediction{SVMpass} == results.minclass & data.testevents == results.minclass) / sum(data.testevents == results.minclass);
        results.(cfg.timecourses{p}).accuracy.chancelevels(SVMpass) = sum(data.testevents == results.minclass)/max(size(data.testevents));
        results.(cfg.timecourses{p}).accuracy.normal(SVMpass)= sum(results.(cfg.timecourses{p}).testprediction{SVMpass} == data.testevents) / length(data.testevents);
        results.(cfg.timecourses{p}).actualuseTest{SVMpass}  = data.actualuseTest;
        results.(cfg.timecourses{p}).actualuseTrain{SVMpass} = data.actualuseTrain;
        results.(cfg.timecourses{p}).TestrunsIndex{SVMpass}  = cfg.TestrunsIndex{SVMpass};
        results.(cfg.timecourses{p}).usedTimecourse{SVMpass}  = data.(cfg.timecourses{p});
        %results.(cfg.timecourses{p}).weights{SVMpass}          = results.
        
        % CALCULATE results.nodetected = results.nopresent ACCURACY:
        if isfield(results,'scaledtestprediction')
            results.(cfg.timecourses{p}).accuracy.normalized(SVMpass) = sum(results.(cfg.timecourses{p}).scaledtestprediction{SVMpass} == 1 & data.testevents == 1) / sum(data.testevents == 1);
            additionaltext = [', raw accuracy: ' mean(results.(cfg.timecourses{p}).accuracy.raw) ', normalized: ' mean(results.(cfg.timecourses{p}).accuracy.normalized) ', chance: ' results.(cfg.timecourses{p}).accuracy.chancelevels];
        else additionaltext = [];
        end
    end
    
    % Print some preliminairy results
    disp([cfg.subject ', ' data.ROI ', timecourse: ' cfg.timecourses{p} ' class ' num2str(results.minclass) ' actual events: ' num2str(mean(results.(cfg.timecourses{p}).nodetected)) ', predicted: ' num2str(mean(results.(cfg.timecourses{p}).nopresent)) ', normal accuracy: ' num2str(mean(results.(cfg.timecourses{p}).accuracy.normal)) additionaltext]);
    if cfg.plotSVMresults
        if ishandle(15); close(15); end;
        pause(0.5);
        SVMplotSVM(cfg,results,cfg.timecourses{p})
        [g gg ggg]=mkdir([cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM results/figures ' cfg.timestarted]);
        disp('Saving plot to file')
        tic;
        saveas(15,[cfg.datadir cfg.experiment{1} '/' cfg.subject '/SVM results/figures ' cfg.timestarted '/' cfg.infotext '.jpg'])
        toc;
    end
end
fclose('all');
end

function matrix = bandpass(cfg,matrix)
%% Applies bandpass
dim = size(matrix);

f = (1/(cfg.Tr/1000))/2*linspace(0,1,ceil(dim(2)/2+0.5));
L1= dim(2)/2+0.5*not(mod(dim(2),2))+0.5;
L2=ceil(dim(2)/2+not(mod(dim(2),2)*0.5))+mod(dim(2),2);

points(2:L1)             = f(1:L1-1)>=cfg.bandpass(1)&f(1:L1-1)<=cfg.bandpass(2);
points(L2:dim(2))        = fliplr(points(2:end));
points(1)                = 1;

for n=1:dim(1) 
    for m=1:dim(3)
        Y               = fft(squeeze(matrix(n,:,m)));
        matrix(n,:,m)   = ifft(Y.*points);
    end; 
end;
end

function matrix = supersmooth(cfg,matrix)
%% Makes the timecourses silky smooth
dim = size(matrix);

smoothspan=round(cfg.smoothing/(cfg.Tr/1000));

if smoothspan <3
    disp('No smoothing, span too small');
end

for n=1:dim(1) 
    for m=1:dim(3)
        matrix(n,:,m) = smooth_backup(squeeze(matrix(n,:,m)),smoothspan);
    end; 
end;

end

function SVMplotSVM(cfg,results,timecoursevariant)
%% Plot Switches
if cfg.plotSVMresults
    fid = figure(15);
    
    for n = 1 : length(results.(timecoursevariant).testevents)
        subplot('Position',[0,(n-1)/length(results.(timecoursevariant).testevents),1,1/length(results.(timecoursevariant).testevents)]);
        hold on
        tf1 = results.(timecoursevariant).testevents{n}==sign(results.(timecoursevariant).predTimecourse{n});
        tf2 = results.(timecoursevariant).testevents{n}~=sign(results.(timecoursevariant).predTimecourse{n});
        im = zeros(1, length(tf1), 3);
        im(1,tf1,2) = 0.5;
        im(1,tf2,1) = 0.5;
        image(im)
        figpos=get(fid,'position');
        uicontrol('position',[0 0 figpos(3) 12],'style','text','string',cfg.infotext,'HorizontalAlignment','left')
        
        plot(results.(timecoursevariant).predTimecourse{n}/max(abs(results.(timecoursevariant).predTimecourse{n}))*0.5+1,'color',[1 1 1])
        line([1 length(results.(timecoursevariant).testevents{n})],[1 1],'color',[1 1 1]);
        axis tight
        axis off
        hold off
    end
    
end

end

function output=SVMjoincell(somecell, cellstojoin)

%% Joins cells and assigns them to test or train

%input: somecell = n cells containing time x voxelvalue?

%output: matrix with n columns, containing the cells of somecell as rows,
%        but only these cells that... WTF DOES THIS CODE DO?!?


% Makes sure no mistakes can be made selecting events and matrices
output                  = [];

if length(somecell)~=size(cellstojoin,2)
    error('Cells to join and somecell have to be equal size');
end

output=[];
for n=1:length(somecell)
    for m = 1:size(somecell{n},2)
        if cellstojoin(m,n)==1
            output(1:size(somecell{n},1),end+1) = somecell{n}(:,m);
        end
    end
end

end

%% caclulating beta in one single regresssion
function betas = findBetas(matrix,timecourse,cfg)
    p = [6 16 1 1 6 0 32];      % Is default
    hrf = spm_hrf_backup(cfg.Tr/1000, p);

    hrfcourse1=conv(double(timecourse== 1),hrf);
    hrfcourse1=hrfcourse1(1:length(timecourse));    
    hrfcourse2=conv(double(timecourse==-1),hrf);
    hrfcourse2=hrfcourse2(1:length(timecourse));

    betas = zeros(2,size(matrix,1));
    for n = 1:size(matrix,1)
        actualusestart = find(timecourse~=timecourse(1),1);
        %discard all points before first switch
        p = regress(matrix(n,actualusestart:end)',[ones(size(hrfcourse1(actualusestart:end))); hrfcourse1(actualusestart:end); hrfcourse2(actualusestart:end)]');
        betas(:,n) = p(2:3);
        
        if (n==10 && cfg.plottimecourses == 1)
           figure(12);
           subplot(2,1,1);
           plot([matrix(n,:); hrfcourse1; hrfcourse2]');
           legend('BOLD','HRFpred1','HRFpred2','Location','EastOutside','Orientation','vertical');
           
           subplot(2,1,2);
           plot([matrix(n,:); (p(1)+ p(2)*hrfcourse1 + p(3)*hrfcourse2)]');
           legend('BOLD','regression','Location','EastOutside','Orientation','vertical');
           
           title(['regression = ' num2str(p(1)) ' + ' num2str(p(2)) ' * HRFpred1 + ' num2str(p(3)) ' * HRFpred2']);
        end
    end
end

%% calculating beta with separate regressions
function betas = findBetas2(matrix,timecourse,cfg)
    p = [6 16 1 1 6 0 32];      % Is default
    hrf = spm_hrf_backup(cfg.Tr/1000, p);

    hrfcourse1=conv(double(timecourse== 1),hrf);
    hrfcourse1=hrfcourse1(1:length(timecourse));    
    hrfcourse2=conv(double(timecourse==-1),hrf);
    hrfcourse2=hrfcourse2(1:length(timecourse));

    betas = zeros(2,size(matrix,1));
    for n = 1:size(matrix,1)
        actualusestart = find(timecourse~=timecourse(1),1);
        %discard all points before first switch
        p1 = regress(matrix(n,actualusestart:end)',[ones(size(hrfcourse1(actualusestart:end))); hrfcourse1(actualusestart:end)]');
        p2 = regress(matrix(n,actualusestart:end)',[ones(size(hrfcourse1(actualusestart:end))); hrfcourse2(actualusestart:end)]');
        
        
        %p1 = regress(matrix(n,:)',[ones(size(hrfcourse1)); hrfcourse1]');
        %p2 = regress(matrix(n,:)',[ones(size(hrfcourse1)); hrfcourse2]');
        betas(:,n) = [p1(2) p2(2)];
        
        if (n==10 && cfg.plottimecourses == 1)
           figure(12);
           subplot(2,1,1);
           plot([matrix(n,:); hrfcourse1; hrfcourse2]');
           legend('BOLD','HRFpred1','HRFpred2','Location','EastOutside','Orientation','vertical');
           
           subplot(2,1,2);
           plot([matrix(n,:); (p1(1)+ p1(2)*hrfcourse1); (p2(1)+ p2(2)*hrfcourse2)]');
           legend('BOLD','regression1','regression2','Location','EastOutside','Orientation','vertical');
           
           %title(['regression = ' num2str(p(1)) ' + ' num2str(p(2)) ' * HRFpred1 + ' num2str(p(3)) ' * HRFpred2']);
        end
    end
end

function betas = findBetasAVG(matrix,timecourse,cfg)
    betas = zeros(2,size(matrix,1));
    for n = 1:size(matrix,1)
        actualusestart = find(timecourse~=timecourse(1),1);
        actualuselogical = zeros(size(timecourse));
        actualuselogical(actualusestart:end) = 1;
        
        %discard all points before first switch
        %p1 = regress(matrix(n,actualusestart:end)',[ones(size(hrfcourse1(actualusestart:end))); hrfcourse1(actualusestart:end)]');
        %p2 = regress(matrix(n,actualusestart:end)',[ones(size(hrfcourse1(actualusestart:end))); hrfcourse2(actualusestart:end)]');
        
        p1 = mean(matrix(n,actualuselogical & timecourse==1));
        p2 = mean(matrix(n,actualuselogical & timecourse==-1));
        
        %p1 = regress(matrix(n,:)',[ones(size(hrfcourse1)); hrfcourse1]');
        %p2 = regress(matrix(n,:)',[ones(size(hrfcourse1)); hrfcourse2]');
        betas(:,n) = [p1 p2];
        
%         if (n==10 && cfg.plottimecourses == 1)
%            figure(12);
%            subplot(2,1,1);
%            plot([matrix(n,:); hrfcourse1; hrfcourse2]');
%            legend('BOLD','HRFpred1','HRFpred2','Location','EastOutside','Orientation','vertical');
%            
%            subplot(2,1,2);
%            plot([matrix(n,:); (p1(1)+ p1(2)*hrfcourse1); (p2(1)+ p2(2)*hrfcourse2)]');
%            legend('BOLD','regression1','regression2','Location','EastOutside','Orientation','vertical');
%            
%            %title(['regression = ' num2str(p(1)) ' + ' num2str(p(2)) ' * HRFpred1 + ' num2str(p(3)) ' * HRFpred2']);
%         end
    end
end